/*                SPSS-Macro TetCorr                     */
/*                                                       */
/*          (Version 2.3; D. Enzmann, 2007)              */
/*                                                       */
/* TetCorr computes the tetrachoric correlation of two   */
/* dichotomous variables and its significance based on a */
/* modified macro written by by Jeffrey S. Kane          */
/* (Associate Professor, Chinese University of Hong Kong */
/* Kong, October20, 1998) that was translated and        */
/* modified from a FORTRAN program in Applied Statistics */
/* (1977), Vol 26, No 3.                                 */
/*                                                       */
/* The variables need not to be coded as 0/1, but the    */
/* codes for 'yes/true' and 'no/false' have to be the    */
/* same for both variables, otherwise TetCorr will       */
/* produce meaningless results.                          */
/*                                                       */
/* If the number of cases is more than 12,0000 you have  */
/* to increase the default value of MXLOOPS to the       */
/* the number of cases in your current working file      */
/* after you installed the macro via INCLUDE file (if    */
/* not, SPSS will most likely issue an error  message).  */
/*                                                       */
/* Version 2.1 did compute wrong estimates if one cell   */
/* frequency was zero. This has been corrected in        */
/* version 2.2, additionally a warning will be issued.   */
/*                                                       */
/* Example:                                              */
/* Two variables x and y are coded '1' for 'yes/true'    */
/* and '2' for 'no/false'. Your data file consists of    */
/* about 22800 cases. First, use the SET MXLOOP command  */
/* once to increase the default value of 12000 to at     */
/* least 22800:                                          */
/* SET MXLOOP=23000.                                     */
/* Next, you can call TetCorr with x and y as parameters */
/* TETCORR x y .                                         */
/*                                                       */
/* A full example can be found in the file EXAMP_R.SPS . */
/*                                                       */
/* For a more exact estimation of the tetrachoric        */
/* correlation in case low cell counts (< 5) or extreme  */
/* asymetric marginal counts see Divgy (1979) and Kirk   */
/* (1973). Tables for estimating the significance in     */
/* case of extreme asymetric marginal counts more        */
/* exactly can be found in Zalinski et al. (1979).       */
/*                                                       */
/* References:                                           */
/* Divgi, D.R. (1979). Calculation of the tetrachoric    */
/*    correlation coefficient. Psychometrika, 44,        */
/*    169-172.                                           */
/* Kirk, D.B. (1973). On the numerical approximation of  */
/*    the bivariate normal (tetrachoric) correlation     */
/*    coefficient. Psychometrika, 38, 259-268.           */
/* Zalinski, J., Abrahams, N.M., & Alf, E.jr. (1979).    */
/*    Computing tables for the tetrachoric correlation   */
/*    coefficient. Educational and Psychological         */
/*    Measurement, 39, 267-275.                          */
/* ----------------------------------------------------- */.

set printback on
   /mxloop=120000.
PRESERVE.
set printback off.

define tetcorr ( !positional !tokens(1)
                /!positional !tokens(1))
matrix.
get raw
   /file=*
   /variables=!1 !2
   /names = vname
   /missing omit.
compute min=mmin(raw).
compute max=mmax(raw).
compute cells=make(1,4,0).
loop #i=1 to nrow(raw).
do if (raw(#i,1)=min and raw(#i,2)=min).
- compute cells(1,1)=cells(1,1)+1.
else if (raw(#i,1)=min and raw(#i,2)=max).
- compute cells(1,2)=cells(1,2)+1.
else if (raw(#i,1)=max and raw(#i,2)=min).
- compute cells(1,3)=cells(1,3)+1.
else if (raw(#i,1)=max and raw(#i,2)=max).
- compute cells(1,4)=cells(1,4)+1.
end if.
end loop.

compute X = {0.9972638618;
             0.9856115115;
             0.9647622556;
             0.9349060759;
             0.8963211558;
             0.8493676137;
             0.7944837960;
             0.7321821187;
             0.6630442669;
             0.5877157572;
             0.5068999089;
             0.4213512761;
             0.3318686023;
             0.2392873623;
             0.1444719616;
             0.0483076657}.
compute W = {0.0070186100;
             0.0162743947;
             0.0253920653;
             0.0342738629;
             0.0428358980;
             0.0509980593;
             0.0586840935;
             0.0658222228;
             0.0723457941;
             0.0781938958;
             0.0833119242;
             0.0876520930;
             0.0911738787;
             0.0938443991;
             0.0956387201;
             0.0965400885}.

compute rlimit= 0.9999.
compute rcut  = 0.95.
compute uplim = 5.0.
compute const = 1E-36.
compute chalf = 1E-18.
compute conv  = 1E-8.
compute citer = 1E-10.
compute niter = 500.
* compute niter = 600000.

/* STRING errmsg1 (A55) / errmsg2 (A60). */
/* string #var1 (A8) / #var2 (A8).       */
/* compute #var1 = !quote(!1).           */
/* compute #var2 = !quote(!2).           */

   /* Initialize variables. */

Compute RTETRA = 0.
Compute SDZERO = 0.
Compute SDR    = 0.
Compute ITYPE  = 0.
Compute endrun = 0.

   /*Check if any cell frequency is negative. */
   /* Proceed if no cell frequency is negative. */

do if cells(1,1)>=0 and cells(1,2)>=0 and cells(1,3)>=0 and cells(1,4)>=0.

   /* Neg freq. end. */
   /* Check if any marginal frequency is zero. */

- compute kdelta = 1.
- compute delta  = 0.
- do if (cells(1,1)=0 or cells(1,4)=0).
-   compute kdelta = 2.
- end if.
- do if (cells(1,2)=0 or cells(1,3)=0).
-   compute kdelta = kdelta+2.
- end if.
- do if (kdelta = 4).
-   compute endrun = 2.
- end if.

   /* If KDELTA=4, a row or column frequency is zero. */
   /* Proceed if no row or column frequency is zero. */

- do if (endrun = 0).
   /* zero row/col end. */

-   do if (kdelta = 2).
-     compute delta = 0.5.
-     do if (cells(1,1)=0 and cells(1,4)=0).
-       compute rtetra = -1.
-       compute margtot = rsum(cells).
-     else
-       compute rtetra = -1
-       compute margtot = rsum(cells).
-     end if.
-   end if.

-   do if (kdelta = 3).
-     compute delta = -0.5.
-     do if (cells(1,2)=0 and cells(1,3)=0).
-       compute rtetra = 1.
-       compute margtot = rsum(cells).
-     else.
-       compute rtetra = -1.
-       compute margtot = rsum(cells).
-     end if.
-   end if.

-   do if rtetra <> 0).
-     compute itype = 3.
-   end if.

   /* Store frequencies in ccells. */

-   compute ccells=make(1,4,0).
-   compute ccells(1,1)=cells(1,1) + delta.
-   compute ccells(1,2)=cells(1,2) - delta.
-   compute ccells(1,3)=cells(1,3) - delta.
-   compute ccells(1,4)=cells(1,4) + delta.
-   compute tot = rsum(ccells).

   /* Check if correlation is negative. */

-   DO IF ((CCELLS(1,1)*CCELLS(1,4))-(CCELLS(1,2)*CCELLS(1,3)))=0.
-     compute itype = 4.
-   end if.

   /* Compute probabilities of quadrants and of marginals. */
   /* PROBAA AND PROBAC chosen so that correlation is positive. */
   /* SIGN indicates whether quadrants have been switched. */

-   do if ((ccells(1,1)*ccells(1,4))-(ccells(1,2)*ccells(1,3))) >= 0.
-     compute probaa = ccells(1,1)/tot.
-     compute probac = (ccells(1,1)+ccells(1,3))/tot.
-     Compute sign = 1.
-   END IF.

-   do if ((ccells(1,1)*ccells(1,4))-(ccells(1,2)*ccells(1,3))) < 0.
-     compute probaa = ccells(1,2)/tot.
-     compute probac = (ccells(1,2)+ccells(1,4))/tot.
-     compute sign = -1.
-   end if.

-   compute probab = (ccells(1,1)+ccells(1,2))/tot.

   /* Compute normal deviates for the marginal frequencies. */

/* zac = idf.normal(probac,0,1) */

/* ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3        */
/*                                                              */
/* Produces the normal deviate Z corresponding to a given lower */
/* tail area of P; Z is accurate to about 1 part in 10**16.     */.

compute SPLIT1 = 0.425.
compute SPLIT2 = 5.
compute CONST1 = 0.180625.
compute CONST2 = 1.6.

/* Coefficients for p close to 0.5  */

compute A = {3.3871328727963666080,
             1.3314166789178437745E+2,
             1.9715909503065514427E+3,
             1.3731693765509461125E+4,
             4.5921953931549871457E+4,
             6.7265770927008700853E+4,
             3.3430575583588128105E+4,
             2.5090809287301226727E+3}.
compute B = {1.0,
             4.2313330701600911252E+1,
             6.8718700749205790830E+2,
             5.3941960214247511077E+3,
             2.1213794301586595867E+4,
             3.9307895800092710610E+4,
             2.8729085735721942674E+4,
             5.2264952788528545610E+3}.

/* Coefficients for p not close to 0, 0.5 or 1 */

compute C = {1.42343711074968357734,
             4.63033784615654529590,
             5.76949722146069140550,
             3.64784832476320460504,
             1.27045825245236838258,
             2.41780725177450611770E-1,
             2.27238449892691845833E-2,
             7.74545014278341407640E-4}.
compute D = {1.0,
             2.05319162663775882187,
             1.67638483018380384940,
             6.89767334985100004550E-1,
             1.48103976427480074590E-1,
             1.51986665636164571966E-2,
             5.47593808499534494600E-4,
             1.05075007164441684324E-9}.

/* Coefficients for p near 0 or 1 */

compute E = {6.65790464350110377720,
             5.46378491116411436990,
             1.78482653991729133580,
             2.96560571828504891230E-1,
             2.65321895265761230930E-2,
             1.24266094738807843860E-3,
             2.71155556874348757815E-5,
             2.01033439929228813265E-7}.
compute F = {1.0,
             5.99832206555887937690E-1,
             1.36929880922735805310E-1,
             1.48753612908506148525E-2,
             7.86869131145613259100E-4,
             1.84631831751005468180E-5,
             1.42151175831644588870E-7,
             2.04426310338993978564E-15}.

compute zac=-7.9414444931916800.
do if probac <= 10E-16.
- compute probac = 10E-16.
else if 1-probac <= 10E-16.
- compute probac = 1 - 10E-16.
end if.
compute exitf=0.

compute ifault = 0.
compute q = probac - 0.5.
do If (Abs(q) <= SPLIT1).
- compute r = CONST1-q*q.
* compute zac = q*(((((((A7*r+A6)*r+A5)*r+A4)*r+A3)*r+A2)*r+A1)*r+A0).
* compute zac = zac/(((((((B7*r+B6)*r+B5)*r+B4)*r+B3)*r+B2)*r+B1)*r+1).
- compute t1=A(8).
- compute t2=B(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+A(-#i).
-   compute t2=r*t2+B(-#i).
- end loop.
- compute zac=q*t1/t2.
/* Exit Function */
- compute exitf=1.
end if.
do If (q < 0) and not exitf.
- compute r = probac.
Else if not exitf.
- compute r = 1 - probac.
End If.
do If (r <= 0) and not exitf.
- compute ifault = 1.
- compute zac = 0.
/* Return */
End If.
do if not exitf.
- compute r = sqrt(-ln(r)).
end if.
do If (r <= SPLIT2) and not exitf.
- compute r = r - CONST2.
* compute zac = (((((((C7*r+c6)*r+c5)*r+c4)*r+c3)*r+c2)*r+c1)*r+C0).
* compute zac = zac/(((((((D7*r+D6)*r+d5)*r+d4)*r+d3)*r+d2)*r+d1)*r+1).
- compute t1=C(8).
- compute t2=D(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+C(-#i).
-   compute t2=r*t2+D(-#i).
- end loop.
- compute zac=t1/t2.
Else if not exitf.
- compute r = r - SPLIT2.
* compute zac = (((((((E7*r+E6)*r+E5)*r+E4)*r+E3)*r+E2)*r+E1)*r+E0).
* compute zac = zac/(((((((F7*r+F6)*r+F5)*r+F4)*r+F3)*r+F2)*r+F1)*r+1).
- compute t1=E(8).
- compute t2=F(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+E(-#i).
-   compute t2=r*t2+F(-#i).
- end loop.
- compute zac=t1/t2.
End If.
do If (q < 0) and not exitf.
- compute zac = -zac.
End If.

/* zab = idf.normal(probab,0,1) */

compute zab=-7.9414444931916800.
do if probab <= 10E-16.
- compute probab = 10E-16.
else if 1-probab <= 10E-16.
- compute probab = 1 - 10E-16.
end if.
compute exitf=0.

compute ifault = 0.
compute q = probab - 0.5.
do If (Abs(q) <= SPLIT1).
- compute r = CONST1-q*q.
* compute zab = q*(((((((A7*r+A6)*r+A5)*r+A4)*r+A3)*r+A2)*r+A1)*r+A0).
* compute zab = zab/(((((((B7*r+B6)*r+B5)*r+B4)*r+B3)*r+B2)*r+B1)*r+1).
- compute t1=A(8).
- compute t2=B(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+A(-#i).
-   compute t2=r*t2+B(-#i).
- end loop.
- compute zab=q*t1/t2.
/* Exit Function */
- compute exitf=1.
end if.
do If (q < 0) and not exitf.
- compute r = probab.
Else if not exitf.
- compute r = 1 - probab.
End If.
do If (r <= 0) and not exitf.
- compute ifault = 1.
- compute zab = 0.
/* Return */
End If.
do if not exitf.
- compute r = sqrt(-ln(r)).
end if.
do If (r <= SPLIT2) and not exitf.
- compute r = r - CONST2.
* compute zab = (((((((C7*r+c6)*r+c5)*r+c4)*r+c3)*r+c2)*r+c1)*r+C0).
* compute zab = zab/(((((((D7*r+D6)*r+d5)*r+d4)*r+d3)*r+d2)*r+d1)*r+1).
- compute t1=C(8).
- compute t2=D(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+C(-#i).
-   compute t2=r*t2+D(-#i).
- end loop.
- compute zab=t1/t2.
Else if not exitf.
- compute r = r - SPLIT2.
* compute zab = (((((((E7*r+E6)*r+E5)*r+E4)*r+E3)*r+E2)*r+E1)*r+E0).
* compute zab = zab/(((((((F7*r+F6)*r+F5)*r+F4)*r+F3)*r+F2)*r+F1)*r+1).
- compute t1=E(8).
- compute t2=F(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+E(-#i).
-   compute t2=r*t2+F(-#i).
- end loop.
- compute zab=t1/t2.
End If.
do If (q < 0) and not exitf.
- compute zab = -zab.
End If.

-   compute ss  = exp(-0.5*(zac**2 + zab**2))/(8*artan(1)).

   /* If RTETRA is 0.0, 1.0, or -1.0, transfer to compute SDZERO. */

-   do if not((rtetra <> 0) or (itype > 0)).

   /* else transfer  to line 85 for r = 0 or 1, -1. */
   /* When marginals are equal, cosine evaluation is used. */

-     do if not(cells(1,1)=cells(1,4) and cells(1,2)=cells(1,3)).

   /* else transfer to line 60. */
   /*Initial estimate of correlation is Yule's Y.*/

-       compute rr = (sqrt(ccells(1,1)*ccells(1,4))
                     -sqrt(ccells(1,2)*ccells(1,3)))**2/
                      abs((ccells(1,1)*ccells(1,4))
                         -(ccells(1,2)*ccells(1,3))).
-       compute iter = 0.

  /* If RR exceeds RCUT, Gaussian quadrature is used. */

-       do if not(rr > rcut).

  /* Tetrachoric series is computed. */
  /* Initialization. */

-         compute endloop = 0.
-         loop.
-           compute term = 1.
-           do if (term = 1).
-             compute va = 1.
-             compute vb = zac.
-             compute wa = 1.
-             compute wb = zab.
-             compute iterm = 0.
-             compute sum1  = probab*probac.
-             compute deriv = 0.
-             compute sr = ss.
-           end if.
-           compute niter2=0.
-           LOOP.
-           compute niter2=niter2+1.
-           DO IF NOT(ABS(SR) > CONST).

  /* Rescale terms to avoid overflows and underflows. */

-             COMPUTE SR = SR/CONST.
-             COMPUTE VA = VA*CHALF.
-             COMPUTE VB = VB*CHALF.
-             COMPUTE WA = WA*CHALF.
-             COMPUTE WB = WB*CHALF.
-           END IF.

  /* Form sum and derivative of series. */

-           Compute DR  = SR*VA*WA.
-           Compute SR  = ((SR*RR)/TERM).
-           Compute COF = SR*VA*WA.

  /* ITERM counts no. of consecutive terms less than CONV. */

-           Compute ITERM = ITERM + 1.
-           do IF (ABS(COF) > CONV).
-             compute ITERM = 0.
-           end if.
-           Compute SUM1  = SUM1 + COF.
-           Compute DERIV = DERIV + DR.
-           Compute VAA   = VA.
-           Compute WAA   = WA.
-           Compute VA    = VB.
-           Compute WA    = WB.
-           Compute VB    = (ZAC*VA) - (TERM*VAA).
-           Compute WB    = (ZAB*WA) - (TERM*WAA).
-           Compute TERM  = TERM + 1.
-         END LOOP IF NOT((ITERM < 2) OR (TERM < 6)).
*                      OR (Niter2 > niter).

  /*Check if iteration converged. */

* -         print niter2 /formats F13.0.

-         DO IF (ABS(SUM1 - PROBAA) > CITER)
                or (niter2 > niter).

  /* If too many iterations, run is terminated. */

-           Compute ITER = ITER + 1.
-           DO IF (ITER < NITER).
*                   and not (Niter2 > niter).
-             Compute DELTA  = (SUM1-PROBAA)/DERIV.
-             Compute RRPREV = RR.
-             Compute RR     = RR-DELTA.
-             do IF (ITER=1).
-               compute RR = RR + 0.5 + DELTA.
-             end if.
-             do IF (RR > RLIMIT).
-               compute RR = RLIMIT.
-             end if.
-             do IF (RR < 0).
-               compute RR = 0.
-             end if.
-           ELSE.
-             Compute endrun  = 3.
-             Compute ENDLOOP = 1.
-           END IF.
-           ELSE.
-             Compute ITYPE   = TERM.
-             Compute ENDLOOP = 1.
-           END IF.
-         END LOOP IF (ENDLOOP > 0).
-       ELSE .

  /* line 40 routine. */
  /* Gaussian Quadrature. */
  /* Initialization, if this is first iteration. */

-         DO IF NOT(ITER > 0).
-           Compute SUM1   = PROBAB*PROBAC.
-           Compute RRPREV = 0.
-         END IF.
-         Compute SUMPRV = PROBAB-SUM1.
-         Compute PROB   = BB/TOT.
-         do IF (SIGN = -1).
-           compute PROB = AA/TOT.
-         end if.
-         Compute ITYPE = 1.

  /* Loop to find estimate of correlation. */
  /* Computation of integral (sum) by quadrature. */

-         LOOP.
-           Compute RRSQ = SQRT(1 - RR**2).
-           Compute AMID = 0.5*(UPLIM+ZAC).
-           Compute XLEN = UPLIM - AMID.
-           Compute SUM1   = 0.
-           Compute endrun = 0.
-           Compute LOOPEND= 0.
-           LOOP #j = 1 TO 16.
-             Compute XLA = AMID + (X(#j,1)*XLEN).
-             Compute XLB = AMID - (X(#j,1)*XLEN).

  /* To avoid underflows, TEMPA and TEMPB are used. */

-             Compute TEMPA = (ZAB-(RR*XLA))/RRSQ.
-             do IF (TEMPA >= -6).
-               compute SUM1 =
                  SUM1 + W(#j,1)*EXP(-0.5*XLA**2)*CDFNORM(TEMPA).
-             end if.
-             Compute TEMPB = (ZAB-(RR*XLB))/RRSQ.
-             do IF (TEMPB >= -6).
-               compute SUM1 =
                  SUM1 + W(#j,1)*EXP(-0.5*XLB**2)*CDFNORM(TEMPB).
-             end if.
-           END LOOP.
-           Compute SUM1 = (SUM1*XLEN)/sqrt(8*artan(1)).

  /* Check if iteration has converged. */

-           DO IF NOT(ABS(PROB-SUM1) <= CITER).
-             Compute ITER = ITER + 1.

  /* If too many iterations, run is terminated. */

-             do IF (ITER >= NITER).
-               compute endrun = 3.
-             end if.
-             DO IF (endrun = 0).
-               Compute RREST = ((PROB-SUM1)*RRPREV-(PROB-SUMPRV)*RR)/
                                 (SUMPRV-SUM1).

  /* Is estimate positive and less than upper limit?. */

-               do IF (RREST > RLIMIT).
-                 compute RREST = RLIMIT.
-               end if.
-               do IF (RREST < 0).
-                 compute RREST = 0.
-               end if.
-               Compute RRPREV = RR.
-               Compute RR     = RREST.
-               Compute SUMPRV = SUM1.

  /* If estimate has same value on two iterations, stop iteration. */

-             END IF.
-           ELSE.
-             Compute LOOPEND = 1.
-           END IF.
-         END LOOP IF ((RR = RRPREV) OR (LOOPEND > 0) OR (endrun > 0)).
-       END IF.
-     ELSE .

  /*Line 60 routine for when all marginals are equal. */

-       Compute RR    = -COS(8*artan(1)*PROBAA).
-       Compute ITYPE = 2.
-     END IF.
-     DO IF ((endrun = 0) and (abs(rtetra) <> 1)).
-       Compute RTETRA = RR.

  /*Compute SDR using McNemar's (1969) approximation formula 12-7. */

-       Compute MargTot = rsum(cells).
-       Compute P1 = (cells(1,1) + cells(1,3))/MargTot.
-       Compute Q1 = (cells(1,2) + cells(1,4))/MargTot.
-       Compute P2 = (cells(1,1) + cells(1,2))/MargTot.
-       Compute Q2 = (cells(1,3) + cells(1,4))/MargTot.
-       Compute HFactor = 1/sqrt(8*artan(1)).
-       DO IF (P1 <= Q1).
-         Compute H1Z = P1.
-       ELSE.
-         Compute H1Z = Q1.
-       END IF.
-       DO IF (P2 <= Q2).
-         Compute H2Z = P2.
-       ELSE.
-         Compute H2Z = Q2.
-       END IF.
-       Compute H1 = HFactor*(EXP(-H1Z/2)).
-       Compute H2 = HFactor*(EXP(-H2Z/2)).
-       Compute STDERRT = (SQRT(P1*Q1*P2*Q2)/(H1*H2*SQRT(MargTot)))*
                          (SQRT((1 - RTETRA**2)*
                          (1-((ARSIN(RTETRA)/90)**2)))).
-     END IF.
-   ELSE.

   /* Line 85 destination for RTETRA = 1.0, 0.0, or -1.0. */

-     Compute SDZERO = SQRT(((ccells(1,1)+ccells(1,2))*
                             (ccells(1,1)+ccells(1,3))*
                             (ccells(1,2)+ccells(1,4))*
                             (ccells(1,3)+ccells(1,4)))/TOT)/
                             (TOT**2 *SS).
-     DO IF (RTETRA = 0).
-       Compute STDERRT = SDZERO.
-     ELSE.
-       Compute STDERRT = -9999.
-     END IF.
-   END IF.
- ELSE .

   /*From zero row / column ender. */

-   Compute endrun = 2.
- END IF.
ELSE .

   /*From negative frequency ender. */

-   Compute endrun = 1.
END IF.

 /* Line 92 error message routine. */

DO IF (endrun > 0).
- compute errmsg = {' '}.
- compute four= make(2,2,0).
- compute four(1,1) = cells(1,1).
- compute four(1,2) = cells(1,2).
- compute four(2,1) = cells(1,3).
- compute four(2,2) = cells(1,4).
- print {'    x = ',vname(1),'    y = ',vname(2)}
       /title 'Variables:'
       /format A8.
- print four
       /format = F8.0
       /title = 'cell counts:'
       /rlabels = 'x-low  ','x-high '
       /clabels = 'y-low  ','y-high '.
- print {errmsg}
       /format A1
       /title 'Error: The tetrachoric correlation could not be computed.'.
- do IF (endrun = 1).
-   print {errmsg}
         /format A1
         /title '       Reason: One or more cell frequencies were negative.'.
- end if.
- do IF (endrun = 2).
-   print {errmsg}
         /format A1
         /title '       Reason: The table has a zero row or zero column.'.
- end if.
- do IF (endrun = 3).
-   print {errmsg}
         /format A1
         /title '       Reason: Result did not converge within max iterations.'.
- end if.
ELSE.
- do if (abs(rtetra) <> 1).
-   compute rtetra=sign*abs(rtetra).
- end if.
- compute px=(cells(1,1)+cells(1,2))/MargTot.
- compute py=(cells(1,1)+cells(1,3))/MargTot.
- compute qx=(1-px).
- compute qy=(1-py).

/* z1 = idf.normal(px,0,1) */

compute z1=-7.9414444931916800.
do if px <= 10E-16.
- compute px = 10E-16.
else if 1-px <= 10E-16.
- compute px = 1 - 10E-16.
end if.
compute exitf=0.

compute ifault = 0.
compute q = px - 0.5.
do If (Abs(q) <= SPLIT1).
- compute r = CONST1-q*q.
* compute z1 = q*(((((((A7*r+A6)*r+A5)*r+A4)*r+A3)*r+A2)*r+A1)*r+A0).
* compute z1 = z1/(((((((B7*r+B6)*r+B5)*r+B4)*r+B3)*r+B2)*r+B1)*r+1).
- compute t1=A(8).
- compute t2=B(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+A(-#i).
-   compute t2=r*t2+B(-#i).
- end loop.
- compute z1=q*t1/t2.
/* Exit Function */
- compute exitf=1.
end if.
do If (q < 0) and not exitf.
- compute r = px.
Else if not exitf.
- compute r = 1 - px.
End If.
do If (r <= 0) and not exitf.
- compute ifault = 1.
- compute z1 = 0.
/* Return */
End If.
do if not exitf.
- compute r = sqrt(-ln(r)).
end if.
do If (r <= SPLIT2) and not exitf.
- compute r = r - CONST2.
* compute z1 = (((((((C7*r+c6)*r+c5)*r+c4)*r+c3)*r+c2)*r+c1)*r+C0).
* compute z1 = z1/(((((((D7*r+D6)*r+d5)*r+d4)*r+d3)*r+d2)*r+d1)*r+1).
- compute t1=C(8).
- compute t2=D(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+C(-#i).
-   compute t2=r*t2+D(-#i).
- end loop.
- compute z1=t1/t2.
Else if not exitf.
- compute r = r - SPLIT2.
* compute z1 = (((((((E7*r+E6)*r+E5)*r+E4)*r+E3)*r+E2)*r+E1)*r+E0).
* compute z1 = z1/(((((((F7*r+F6)*r+F5)*r+F4)*r+F3)*r+F2)*r+F1)*r+1).
- compute t1=E(8).
- compute t2=F(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+E(-#i).
-   compute t2=r*t2+F(-#i).
- end loop.
- compute z1=t1/t2.
End If.
do If (q < 0) and not exitf.
- compute z1 = -z1.
End If.

/* z2 = idf.normal(py,0.1) */

compute z2=-7.9414444931916800.
do if py <= 10E-16.
- compute py = 10E-16.
else if 1-py <= 10E-16.
- compute py = 1 - 10E-16.
end if.
compute exitf=0.

compute ifault = 0.
compute q = py - 0.5.
do If (Abs(q) <= SPLIT1).
- compute r = CONST1-q*q.
* compute z2 = q*(((((((A7*r+A6)*r+A5)*r+A4)*r+A3)*r+A2)*r+A1)*r+A0).
* compute z2 = z2/(((((((B7*r+B6)*r+B5)*r+B4)*r+B3)*r+B2)*r+B1)*r+1).
- compute t1=A(8).
- compute t2=B(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+A(-#i).
-   compute t2=r*t2+B(-#i).
- end loop.
- compute z2=q*t1/t2.
/* Exit Function */
- compute exitf=1.
end if.
do If (q < 0) and not exitf.
- compute r = py.
Else if not exitf.
- compute r = 1 - py.
End If.
do If (r <= 0) and not exitf.
- compute ifault = 1.
- compute z2 = 0.
/* Return */
End If.
do if not exitf.
- compute r = sqrt(-ln(r)).
end if.
do If (r <= SPLIT2) and not exitf.
- compute r = r - CONST2.
* compute z2 = (((((((C7*r+c6)*r+c5)*r+c4)*r+c3)*r+c2)*r+c1)*r+C0).
* compute z2 = z2/(((((((D7*r+D6)*r+d5)*r+d4)*r+d3)*r+d2)*r+d1)*r+1).
- compute t1=C(8).
- compute t2=D(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+C(-#i).
-   compute t2=r*t2+D(-#i).
- end loop.
- compute z2=t1/t2.
Else if not exitf.
- compute r = r - SPLIT2.
* compute z2 = (((((((E7*r+E6)*r+E5)*r+E4)*r+E3)*r+E2)*r+E1)*r+E0).
* compute z2 = z2/(((((((F7*r+F6)*r+F5)*r+F4)*r+F3)*r+F2)*r+F1)*r+1).
- compute t1=E(8).
- compute t2=F(8).
- loop #i=-7 to -1.
-   compute t1=r*t1+E(-#i).
-   compute t2=r*t2+F(-#i).
- end loop.
- compute z2=t1/t2.
End If.
do If (q < 0) and not exitf.
- compute z2 = -z2.
End If.

- COMPUTE dx = EXP(-0.5*z1**2)/sqrt(8*artan(1)).
- COMPUTE dy = EXP(-0.5*z2**2)/sqrt(8*artan(1)).
- compute sigmat = sqrt((px*qx*py*qy)/MargTot)/(dx*dy).
- compute Z = rtetra/sigmat.
- compute sig = (1-cdfnorm(abs(z)))*2.
- compute four= make(2,2,0).
- compute four(1,1) = cells(1,1).
- compute four(1,2) = cells(1,2).
- compute four(2,1) = cells(1,3).
- compute four(2,2) = cells(1,4).
- print {'    x = ',vname(1),'    y = ',vname(2)}
       /title 'Variables:'
       /format A8.
- print four
       /format = F8.0
       /title = 'cell counts:'
       /rlabels = 'x-low  ','x-high '
       /clabels = 'y-low  ','y-high '.
- do if (abs(rtetra) <> 1).
-   print {margtot,rtetra,sig}
         /format F13.5
         /title 'tetrachoric correlation:'
         /clabels = 'N','r','p(2-sided)'.
- else
-   print {margtot,rtetra}
         /format F13.5
         /title 'tetrachoric correlation:'
         /clabels = 'N','r'.
-   compute errmsg = {' '}.
-   print {errmsg}
         /format A1
         /title 'WARNING: Zero cell counts; p not computed!'.
- end if.
END IF.
END MATRIX.

!enddefine.
restore.

/* ---------------------------------------------- */.
/* TetCorr is called by:                          */.
/*                                                */.
/* TETCORR var1 var2.                             */.
/*                                                */.
/* Remember to set the MXLOOP setting to the      */.
/* number of cases in your data file the first    */.
/* time you call TETCORR by using                 */.
/* SET MXLOOP=nnn.                                */.
/* (with nnn >= number of cases)                  */.
/* ---------------------------------------------- */.
